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Abstract 

We demonstrate a method which allows the stochastic modelling of quantum 
systems for which the generalised Fokker-Planck equation in the phase space 
contains derivatives of higher than second order. This generalises quantum 
stochastics far beyond the quantum-optical paradigm of three and four-wave 
mixing problems to which these techniques have so far only been applicable. 
To verify our method, we model a full Wigner representation for the optical 
parametric oscillator, a system where the correct results are well known and 
can be obtained by other methods. 
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I. INTRODUCTION 



Coherent-state, or phase-space, path integrals |l] are so far the only variety of path 
integrals that can be evaluated numerically for real time problems. For interactions which 
are no more than quadratic in creation or annihilation operators, the measure over the paths 
can be described in terms of stochastic differential equations (SDEs, Langevin equations) 
for the paths Unlike, e.g., the Feynman path integral, such real-time phase-space 

path integrals are characterised by a constructively defined positive measure, and can be 
calculated by simulating the corresponding Langevin equations. These techniques have been 
most successfully used in quantum optics J2|, but have also been applied to the quantum 
dynamics of condensed bosonic atoms @]. Mathematically, however, the existence of an 
SDE for the paths is restricted to problems were the equation for the corresponding pseudo- 
probability distribution is a genuine Fokker-Planck equation (this is the content of Pawula's 
theorem ||). This imposes the above restriction on the interaction Hamiltonian and thus 
confines the class of problems for which the measure of the path integral may be characterised 
in terms of an SDE to three and four-wave mixing problems. (The two-body collisional 
interaction of bosons has the formal structure of four wave mixing.) 

There exist, however, a wide variety of problems where the generalised Fokker-Planck 
equation is of third or higher order and hence Langevin equations may not be derived. As 
an example, a description of nonlinear quantum optical processes in terms of the Wigner 
distribution results in generalised Fokker-Planck equations with third-order derivatives |§. 
A common procedure is to truncate the Wigner equation at second order, which is equivalent 
to using the semiclassical theory of stochastic electrodynamics 0. This procedure necessarily 
discards the deeper quantum aspects of the problem and gives answers at odds with quantum 
mechanics for several systems 0. There are also situations where even a P-representation 
Fokker-Planck equation must also be written in generalised form |J, and more are likely to 
be investigated in the future. It is therefore of considerable interest to generalise methods 
allowing for a constructive characterisation of the measure of the phase-space path integral 
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beyond the very restrictive quantum-optical paradigm. 

Pawula's theorem would at first sight seem to forbid this generalisation. However, on 
closer inspection we see that the restrictions imposed can be evaded by discretising time. 
Full technical details are presented elsewhere || ; the aim of this letter is to demonstrate the 
feasibility of this generalisation, using as a demonstrative example the Wigner representation 
for the optical parametric oscillator (OPO). We chose this system because it is well known 
and the results can be obtained by other methods. We derive a system of stochastic difference 
equations in a doubled phase-space, which is related to the Wigner representation in the 
same way as the well-known positive-P equations [|Hj] are related to the P-representation. 
By analogy, we shall call this representation the positive-W representation. We show that 
this representation gives the correct results for quadrature relaxation, whereas the truncated 
Wigner representation makes noticeably different predictions. 

It should be stressed that no continuous time limit exists for the positive-W equations. 
Unlike the Wiener process, where sampling noise is independent of the time step for a given 
sample size, here sampling noise diverges in the continuous time limit. This is how Pawula's 
theorem is enforced. In practice, however, one is interested in the sampling noise vs time 
step not for a given sample size, but for a given computational time. From this perspective, 
the difference we find is much less dramatic. Although, as the time step is decreased, the 
computational time grows at a faster rate than for conventional stochastic integration, this 
dependence remains polynomial so that the problem stays in the same class of computational 
complexity. 



II. POSITIVE-W REPRESENTATION 

Degenerate optical parametric oscillation is an optical process in which a nonlinear 
medium inside an optical cavity is pumped with light at one frequency and emits light 
at half that frequency. The free, damping and pumping Hamiltonians for this system may 
be written in their usual forms |TTJ while we choose the nonlinear coupling constant between 
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the light modes, k, to be real so that the interaction Hamiltonian is 



# int - -y 



t 2 fo _ a 2 6+ . (1) 



The annihilation (creation) operators a (at) and b (¥) annihilate (create) photons at the 
lower and higher frequencies respectively. Proceeding via the usual methods, we can map 
the Hamiltonian onto differential equations for the Wigner and positive-P distributions 0. 
The positive-P representation gives a Fokker-Planck equation and can thus be mapped 
straightforwardly onto a set of coupled SDEs using Ito rules. The equation of motion for 
the Wigner distribution, however, has third-order derivatives and thus has no mapping onto 
stochastic differential equations. 

The positive-P representation can alternatively be derived || by postulating a mapping of 
time-normal (Tjv) averages |L2| of the Heisenberg quantum-field operators, st(t), bt(t),a(t), 



b(t), onto classical averages. For example, 



= a!(t)W(t')a(t")b(t"') . (2) 

The upper bar on the RHS of this relation denotes averaging over the statistics of the four 
stochastic c- number fields, a)(t),tf(t),a(t),b(t) (i.e. phase-space path integration). The 
essence of the positive-P representation is in defining these statistics constructively. The 
Heisenberg equations of motion are mapped onto the SDE's for the fields, while the averaging 
over the initial state of the quantum field (denoted as (• • •)) is mapped onto the distribution 
over the initial conditions. 

The Tjy ordering emerges when we generalise the normal ordering of free-field operators 
to the Heisenberg operators. Generalising the symmetric, or Wigner, ordering of free-field 
operators results in a time- Wigner (Tyy) ordering of the Heisenberg field operators |J. Under 
T/v, the "most recent" creation (annihilation) operator becomes the leftmost (rightmost) in 
the product. The TV-ordering acts in a similar way except that it is symmetrised with 
respect to the creation and annihilation operators: 



t w {mm ■ ■ ■ Kt")} = ~ [t w m 1 ) ■ ■ ■ *(0> m 

+ x(t)T w {m ■ ■ ■ Kt")}} ■ (3) 

Here, x, y, ■ ■ ■ ,z stand for field operators (i.e. a or a^), and t should exceed all other time 
arguments in the product so that t > t', ■ ■ ■ , t" . Equation (|3|) is a recurrence relation which 
defines TV f° r the case of different time arguments. The full definition may be found in || . 

Following the example of the positive-P mapping (fj), we postulate a positive- W mapping 
relating TV-ordered operator averages to classical averages of four stochastic c-number fields, 
a(t), o^(t), (3{t), ffl(t), so that, e.g., 

(T w {k\t)b\t')a{f)b{t m )}) = rf(t)p(t>)a(t")f3(t'"). (4) 

We emphasise that this mapping is distinct from that of Eq. (||D by changing the notation 
of the c-number fields. In full detail, the quantum-field-theoretical (QFT) techniques which 
we used in order to characterise the mapping (H) constructively are described elsewhere 



TJj. For the purposes of this letter, we note that these are a straightforward adaptation 



of similar QFT techniques described in Refs. [rjj,[0J. In JH] , Matsubara-style quantum 
dynamics were mapped onto an imaginary-time SDE; in []15|] , Feynman diagram techniques 
were mapped onto a real-time SDE. Here, we consider dynamics on the Schwinger-Keldysh 
C-contour ||16|| . To include the Wigner representation, one needs a generalisation of Wick's 



theorem to the case of symmetric ordering. This generalisation is quite straightforward 



and results in a Keldysh-style diagram series for the TV-ordered averages. As in |1J,|T5 
propagators in this diagram series are then expressed by the retarded Green's function of the 
free Schrodinger equation, and the whole series is restructured so as to make this Green's 



function a propagator in a new series. This yields a Wyld-type series ||17|| , also termed causal 
series . By applying multiple Hubbard-Stratonovich transformations , we eventually 
arrive at a classical stochastic problem for which this series is a solution. This final step of 
the derivation is of independent interest and is discussed in more detail below. 

In the strict mathematical sense, this derivation fails: no continuous-time process exists 
satisfying Eq. (Q). The way around this problem is to allow the mapping (f|) to hold only 



approximately, and consider stochastic processes in discretised time. We then find the 
following set of stochastic difference equations, (dropping the t-dependence for notational 
simplicity) 

Act = (-71a + na ] (3) At + ^r/iAt 1 / 2 + aiAt 1/3 , 
Aa t = (-7ia t + Kaft) At + ^mlAt 1/2 + a\At 1/3 , 
A(3=[e- l2 f3 - |a 2 ) At + v/^At 1 / 2 + a 2 At^, 

A/3* = (e* - 72 /?t - |«t 2) At + ^* At i/2 + a t At i/3 ; (5) 

where At is the step of time discretisation, Aa(t) = a(t + At) — a(t) (likewise for the other 
field variables) and 771,2(0 ar e independent complex standardised Gaussian noises such that 



m(t) = m(t) = 0, m (t)vi(t') = mi^it') = o, V i(t)vt(t') = V 2(t)vm = <W. (6) 

It should be noted that the 77's are 5 a > (Kronecker) correlated, not 8(t — t') (Dirac) correlated 
and that we have explicitly taken care of the proper powers of At. The 7^, j = 1, 2, are the 
cavity loss rates at each frequency and e represents the classical pump. For the cr's we have 



where £1, £1,^2, £2 are independent complex standardised Gaussian noises (with the same 
properties as 771,772). The other parameters obey the relations 



, rs^ = r^s = 1. (8) 



Within these constraints they can be chosen at will and may in fact even be field and/or 
time-dependent. This freedom can be used to control sampling noise in simulations. 
Comparing equations (H) to the partial differential equation for the ^-function, 



dW(a,P,t) 
di 



d d 
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W(a,P,t), (9) 



we see that there is an obvious one-to-one correspondence between n-th order derivatives 
in (|9|) and terms oc At 1 /™ on the RHS's of Eqs. (pj). More precisely stated, each term in 
(0) corresponds to a particular cumulant |0| of increments in (El). This means that, for ex- 



ample, (71a — na*(3) W specifies Aa = (—71a + na*f3) At, resulting in the contribution 
(—71a + na* (3) At to Aa. Note that the double upper bar is used as a notation to signify 



cumulants. We can now also see that 7i Q^ a * W specifies AaAa^ = 71 At and thus yields 
the contributions 771 1/71 At to Aa and T)ly/jiA~i to Aa^. Finally, the third order derivatives 
are represented by the cr's, with the latter defined in such a way as to have only two nonzero 
third-order cumulants of the increments, 

nAt 



Aa 2 A/5t = Aat 2 A/3 = — . (10) 

The one-to-one correspondence between (pj) and @ suggests that rules may be devised 
for finding coefficients in the positive-W equations, starting from the equation for the W- 
function. (This would however not constitute their derivation: whereas the H^-function 
applies only to same-time symmetrically ordered operator averages, (|J) and (0) cover a 
much wider class of multi-time, time-Wigner ordered averages.) QFT methods may then 
be of much assistance when factorising "noise tensors" such as fllOl). Consider, for example, 
the way in which the cr's were found from the cumulants (^) (while the latter were actually 
obtained using the QFT techniques). Comparing Eqs. (|l0l) to Eqs. (|), we see that the 
(same-time) cr's may be specified by postulating the characteristic function: 



$(Ci, Ci f , C2, C 2 f ) = &°l+Cl<n+to4+G°* = e -< C 2 /8-<c 2 T /8_ (n) 



Although real- valued noises certainly cannot exist which satisfy this definition |19| , they are 
easily constructed as complex noises. To this end, consider a complex Hubbard- Stratonovich 
transformation |18[ : (x,y are arbitrary numbers) 

e xy = ffi e xt+yii*-\e = eXi+yi * L y ^ x ^ + y £*\ . (12) 

J IT I 

Here, ^ is a standardised complex Gaussian noise, with the probability density e - '^' 2 /n. The 
formula in square brackets defines a convenient shorthand; using it, we may write: 



d 2 c 2 /8 =^ cU& + ck^e 2 A cu& + cU^Jm + c2T ] i\Jm. (13) 



We have thus recovered the expressions for the a±, a\ pair; the a\, a% pair is derived similarly. 

The a noises successfully mimic genuine (non-classical) third-order noises, given that 
averages involving their complex conjugates never occur. The latter is indeed the case for 
equations @. With the complex conjugates included, we find nonzero cumulants of arbitrary 
order, exactly as expected for non-Gaussian statistics. Note that the necessity of eliminating 
complex-conjugate noises is exactly the reason why introducing various nonclassical noises 
requires a doubling of the phase space; this is no different from the positive-P representation. 

In order for equations ([5]) to match the mapping (^), the values of cumulants mixing 
the increments in ([5]) with their complex-conjugates are irrelevant, but from a practical 
perspective they do in fact turn out to be very relevant as they affect the sampling noise. 
Consequently the sampling errors can be minimised by using the freedom in the definitions 
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of the noises. In the present case minimisation of the quantity |<7i| + 



+ X [ l^l 2 + 



has a noticeably beneficial effect on the numerical integration, with x being a free parameter 
which can be used to redistribute noise between the two modes. Noting that |£| = y/n, we 
find the minimum at 



P = P ] = — - — Tit , s = s ] = x 1/4 , (14) 
4( X vr) 1/6 

which via fl8|) also fixes the values of r's and g's. 

III. NUMERICAL RESULTS 

By numerical experiments we found that, for values of k — jj — 1 and e = 1.5e c , where 
e c = 7 2 /k, a value of x — 0-33 gave the most stable results. We should note here that we 
are working in a very strong-interaction regime because this is where it is easiest to see the 
differences between the positive-P and truncated Wigner results. This regime is however 
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not necessarily unphysical ||2U||. We begin our integration with initial conditions taken as the 



symmetry-broken semiclassical steady state solutions above threshold [11], with a positive. 



The true physical average is however zero, as a may also be negative. Hence, starting with 
this initial condition, we see a decay of the value of the quadrature average (X a ) = a + ofl 
due to quantum tunneling from positive to negative values. As stated above, it is well known 
that the positive-P and truncated Wigner representations give different predictions for this 
tunnelling. 

We have numerically integrated the equations of motion in the three representations, 
using a standard Euler technique. We find excellent agreement between the positive-P and 
positive-W results, as can be seen in Fig. [I], which shows the short time results for the 
quadrature relaxation in the OPO. The positive-W was averaged over 1.2 x 10 7 trajectories, 
with 2.7 x 10 6 for the positive-P, and 1.7 x 10 6 for the truncated Wigner. This was sufficient 
to ensure excellent convergence over the range plotted. Although the positive-W eventually 
falls victim to enormous sampling errors, where it converges, it reproduces the positive-P 
results almost exactly. The truncated Wigner makes noticeably different predictions. 

In summary, we have developed and described a computational method for numerical 
modelling of processes which result in generalised Fokker-Planck equations with third-order 
derivatives. Although we cannot define a continuous limit of our method as a stochastic 
process, this is not operationally important as interesting systems requiring representation 
with third order noises are likely to be treated numerically. We have successfully demon- 
strated our method for the example of quantum tunneling in the OPO, where the neglect 
of third-order terms in the Wigner representation is known to give erroneous results. The 
success with this method gives confidence that the technique may be used to model processes 
where a P-representation may require third-order derivatives. The real importance of our 
method is that it may be used to extend the use of stochastic integration via the phase-space 
representations beyond the field of quantum optics, allowing the deeply quantum aspects of 
a wider range of systems to be investigated. 
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FIG. 1. Predictions of the positive-W (upper solid line), positive-P (dashed line) and truncated 
Wigner (solid line) representations for quadrature relaxation in the OPO. We can see that the 
positive-P and positive-W solutions are almost indistinguishable. 
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